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Abstract 

A techniqu e is presented for adapting computational meshes used in 
the C2 version of the direct, simulation Monte Carlo method. The physical 
ideas underlying the technique are discussed, and adaptation formulas are 
developed for use on solutions generated from an initial mesh. The effect of 
statistical scatter on adaptation is addressed, and results demonstrate the 
ability of this technique to achieve more accurate results wit hout increasing 
necessary computational resources. 

Introduction 

In recent years, statistical particle methods, such as the direct, sim- 
ulation Monte Carlo method (DSMC), have become a popular approach 
for simulating high-speed, rarefied gas flows. 1-4 However, their widespread 
use is often hampered by heavy demands on computational resources, such 
as memory and CPU time. In addition, investigators must often use their 
initial results as a guideline for obtaining more accurate answers by modi- 
fying initial spatial or temporal resolution. This iterative process increases 
the amount of user time, and often the amount of computational resources, 
necessary to compute an accurate solution. 

The purpose of this investigation is to examine algorithms for automat- 
ically adapting grids to initial DSMC* solutions using a number of seemingly 
conflicting criteria for accuracy. With such procedures, we expect to reduce 
the memory, CPU time, and user time necessary to achieve satisfactory lev- 
els of accuracy. The algorithms examined here are compatible with the G2 
version of the DSMC' method. 
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Discussion 

Prom experience with DSMC simulations, users have developed accu- 
racy criteria with respect, to grid resolution. While it is unnecessary to obey 
these rules everywhere at all times, following them provides some confidence 
that the final results are as accurate as the method permits. 5 

One rule maintains that in viscous layers, shock waves, and shear 
layers, where macroscopic properties undergo strong changes, cell spacing 
along these gradients should not, exceed one local mean free path length, A. 
This ensures that properties transported across these gradients through 
interparticle collisions are captured properly. 

Another rule requires that near surfaces, cell spacing normal to the 
body should also be shorter than A. This is because in certain areas, such 
as in the vicinity of adiabatic surfaces, flowfield gradients may be very small, 
and the first rule is not stringent enough to adequately capture flowfield 
physics. 


Adaptation Methodology 

Intuitively, it appears that functions advantageous for grid adaptation 
should have properties pertinent to important aspects of the flowfield as 
well as having properties advantageous to fulfilling the rules cited above. 
Regarding the first criterion, areas where strong gradients occur are often 
characterized by a non negligible “gradient Knudsen number," A'n rt . 5 The 
reciprocal length scale (/«) _1 of this Knudsen number is comprised of the 
local change in some flow variable tv along its direction of strongest change £ , 
divided by the local value of the flow variable itself 
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Gradient length scales appear to he desirable quantities to use as physical 
length scales for grid adaptation. As just noted however, variables chosen 
for rt, such as density, velocity, or temperature, while important in much of 
the flowfield, may not accurately reflect changes near certain surfaces. To 
satisfy the second criterion, we use A as another physical length scale for 
grid adaptation. Presently A is computed using the equilibrium relation- 
ship. Although this formula is not valid where translational nonequilibrium 
exists, it estimates A well enough for initial adaptation purposes at this 
stage of the investigation. 

Using reciprocals of these physical length scab's, we construct functions 
to adapt initial grids from the solutions they generate. To date, the regions 
comprising the flowfield are adapted in one direction only, using either 
rows or columns of cells in the body-fitted mesh. In Fig. 1, these directions 
are denoted by £ and 7/, respectively. Gradients are calculated using values 
contained in cells and their nearest neighbors in the direct ion of adaptation. 
This approximation creates inaccuracies in terms of obtaining true gradients 
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Fig. 1 Schematic representation 
of DSMC domain mesh. 


with respect to a reference frame for the overall flowfield, hut it provides a 
relative measure of them along the direction of adaptation. Although the 
following procedures are discussed in terms of column adaptation, the same 
applies for row adaptation as well. 

An adaptation function F is generated from the absolute reciprocal val- 
ues of these lengths by computing their cumulative sum up through each 
cell in the direction of the grid adaptation. This function increases rnono- 
torncalJy, and has its steepest slopes where the inverse lengths are greatest. 
Typical behavior of F is shown in Fig. 2c. For column adaptation, F is 
a function of 7/, and for a constant number of cells per column, cell spac- 
ing Ar/ is altered through the adaptation process by inverting the function. 
In this manner, the new cell spacing A if is given by equal increments of F, 
and A jf will be minimized where the slope of F is steepest. A typical value 
in a term for altering A?/ within a column of cells in an / x J matrix is 
computed the following way 


1=1 
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As the subscripts indicate, each column of cells has a separate function F 
for adaptation. If F is normalized by its final value, as in Eq. (2), so that 
values lie between zero to unity, a number of such functions, using different, 
variables for a, may be combined to create a composite function G. Addi- 
tionally, the individual terms of this composite function may he weighted 
so that they dominate G where they are strongest. 

T hese ideas were tested for rarefied flow over the front half of a cylinder. 
Fig. 2a shows the initial grid and geometry used in the simulation, which 
took advantage of Howfield symmetry about the stagnation streamline. The 
freest ream Mach, Reynolds, and Knudsen number values were — 3,1, 
ffrr.no — 88. fg and I\n r rxj = 0.052, respectively, based on cylinder radius, 
and the cylinder’s surface temperature was equal to the freest. ream stagna- 




Fig. 2 Flow over cylinder, M co = 3.1, Re r ,oo — 88.6, A n r 0 o — 0.052, &: 
Tw/Too = 2.92: (a) original mesh, (b) density contours [kg/m ], (c) vari- 
ation of F & G for cells nearest centerline. 


tion temperature (r„ = 300 K, = T 0 = 877 K). Density contours from 
the results are shown in Fig. 2b. A composite function G was developed to 
emphasize tight cell spacing near the surface and through the shock wave 
to adequately capture surface fluxes and flowfleld gradients through the 
shock wave. The function adapted cell spacing At/ nominally normal to 
the cylinder’s surface, and was weighted towards A near the surface and 
towards l n , the gradient length scale based on number density n, near the 
outer boundary. The specific form of G became 
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Fig. 3 Meshes adapted to solution depicted in Fig. 2: (a) adapted to 
Eq. (3) without filtering, (b) same, with slight FFT filtering. 


The weighting factor was controlled by the free parameter an exponent 
whose value was set to unity for linear weighting. The variation of G and 
its constituents with distance away from the cylinder near the stagnation 
streamline are shown in Fig. 2c. 

Hg. 3a shows the adapted mesh resulting from application of Eq. (3) 
to the initial results. Comparing Fig. 3a and Fig. 2a, it is apparent that 
mesh lines are now strongly biased towards the location of the shock wave 
and near the cylinder’s surface in the stagnation region. 

Handling of Statistical Fluctuations 

Fig. 3a also reveals a complication in the grid adaptation process. Note 
the “herringbone” pattern resulting from the statistical nature of particle 
methods, which is accentuated in regions of large l n . Cell-averaged flowfield 
properties have a certain level of statistical scatter, or “noise” associated 
with them. In order to shorten the iteration process, we wish to adapt the 
mesh to solutions with small levels of sampling, which magnifies the effects 
of statistical scatter. This problem is further aggravated by adapting to 
fiowfield gradients, obtained by differencing statistically noisy values. 

Efforts to use least-square fits of polynomials and conic sections to the 
adapted grid lines proved ineffective. Such fits shifted tight grid spacing 
away from certain areas where it was physically needed. We also evaluated 
different filters to reduce the noise in the data set and/or the altered mesh. 
They were: centrally weighted averaging, median filtering, 6 and fast fourier 
transform (FFT) filtering/ Of these, FFT filtering has been the most effec- 
tive at reducing the high-frequency noise in our data sets. However, it can 
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produce low-frequency variations in values of data sets or in grid spacing, 
sometimes shifting tight grid spacing away from where it is needed. Nev- 
ertheless, we retained it as a tool for producing acceptable grids. Fig. 3b 
shows the result of applying FFT filtering to the mesh depicted in Fig. 3a. 

The most successful method for reducing statistical scatter encountered 
thus far resulted from observations concerning F and (/. These functions 
adapt to data in only one direction, and for column adaptation each column 
of cells is adapted independently. There is no regard for the fact that final 
values the denominators receive before normalization vary across columns. 
In Fig. 4, the values of F denominators for each column of the simulation 
depicted in Fig. 2 are shown. The first column represents cells closest 
to the flowfield centerline, and the last column represents those nearest 
the exit plane. Such variation depends on mesh geometry, over which F 
denominators may vary widely. When F or G is used where no strong 
gradients exist, the functions tend to adapt to the spiky noise itself. We 
preferred preserving features of the original grid where this occurs. 

A variable coefficient T(/, j \ / a ) was created for F , equal to the ratio of 
the denominator of F for each column to the maximum value of F denom- 
inators across all columns. The unitary complement of this coefficient then 
multiplies the incremental arc length of the original grid itself. 

The coefficient and the modified function F are given by 


t = l 


A 17 (ij) 
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Fig. 4 Variation of F denominators across columns for solution de- 
picted in Fig. 2. 
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and 

F(mJ-l a ) = r(I,jJ a ) F(m,r,l a )+ (i ~ W (5) 


Since F ranges from zero to one, the impact of individual F variations on 
grid spacing are de-emphasized along columns where their summations are 
relatively small. Again, F functions based on different variables may be 
combined to create composite functions, herein denoted by G. 


Results 

In Fig. 5a, we show an initial grid consisting of four separate regions, 
once again for flow over a cylinder. In this example, the freestream con- 
ditions consisted of M ^ = 6.0, Rc r oo = 684, /Tn r oo = 0.013, and 7"^ 
= 63 K. The surface temperature was T UI = 300 K. Since the overall govern- 
ing Knudsen number places this simulation in the near-continuum regime, 
multiple regions enhanced resolution in the vicinity of the bow shock and 
in the viscous region near the surface. Simulation results are qualitatively 
similar to those depicted in Fig. 2b. We fabricated the grid shown in Fig. 5b 
using the following composite function G } based on l n and A 
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By using A as one of the variables for adaptation, we expected to reduce 
in the number of cells violating the rule regarding grid spacing versus A in 
areas of enhanced resolution. Diagnostics in the adaptation code confirm 
this feature for this case. In the initial solution, the fraction of cells violating 
this rule was 0.1% in Region 1 (near the surface), and 24% in Region 3 
(through the shock wave). It should be noted that few of these cells grossly 
violated the rule, however. After repeating the simulation on an adapted 
grid without FFT smoothing for approximately the same number of steady- 
state timesteps, this rule was violated only 0% and 9% in those respective 
regions. 

Up to this point, the simulations presented here have arguably had 
adequate resolution before undergoing grid adaptation. Our final example 
demonstrates that this adaptation process can increase solution accuracy 
for a given grid with a fixed number of cells and regions even when the initial 
solution suffers from inadequate resolution. This case steins from a series 
of simulations originally run in conjunction with a DSMC* investigation of a 
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Fig. 5 Multiple-region mesli for calculating flow over cylinder, Moo 
= 6.0, Re r> o o = 684, A ?i r>00 = 0.013, &: T w /Too = 4.76. Only every fifth 
“tangential” grid line shown for clarity: (a) original mesh, (b) mesh 
adapted to initial solution using Eq. (6), without filtering. 


shock-shock interaction problem. 1 Data from the interaction problem was 
used to model the shear layer that develops beyond the point where the 
shock waves intersect. The simulation was run for three different levels 
of resolution approximately normal to the flow direction. The medium 
resolution grid has four times as many cells in the y-direction as the coarsest 
mesh, and the highest resolution grid has sixteen times as many cells in that 
direction as the coarsest mesh. It is considered that only the final mesh 
had adequate resolution. 1 Fig. 6 depicts velocity magnitude contours from 
these high-resolution results. 



Fig. 6 Shear flow generated using conditions from separate shock-shock 
interaction problem. Velocity magnitude contours [m/sec] . 
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Fig. 7 Meshes used to calculate flowfleld depicted in Fig. 6: (a) original 
mesh, (b) mesh adapted to initial solution using modified version of 
Eq. (6) with slight FFT filtering. 

In Fig. 7a, we show the initial grid used in the coarsest simulation. 
Investigating the effects of different length scales on adaptation in this 
particular case led us to the grid depicted in Fig. 7b. This was generated 
using a composite function G similar to that described in Equation (6), 
using length scales based on /„ and A. The only difference was to replace 
the factor 

f v{™,j) y 

\V(IJ) ) 

with 1/2 because this flowfield did not require the type of weighting desired 
for the bounded flowfields depicted earlier. The length scale based on A was 
included to dampen noisy contributions from l n . 

One run was generated using the grid adapted to the coarsest mesh 
results. Another run was generated using a grid adapted to the medium 
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Fig. 8 Profiles of <pue> for different meshes used to calculate flow- 
field of Fig. 6 at x = 0: (a) unadapted/coarse, (b) adapted/coarse, 
(c) unadapted/medium-resolution, (d) adapted/medium-resolution, 
and (e) unadapted/higli-resolution. 


resolution mesh using the same function G. Contour plots of the adapted 
grid results appeared virtually identical to their respective original cases. 
However, as depicted in Fig. 8, profiles of the energy flux <pue> at a plane 
across x = 0, a location containing some of the steepest gradients, demon- 
strate the relative effectiveness of the grid adaptation technique. The peak 
flux from the solution on the coarse adapted grid is 15% higher than that 
from the plain coarse grid, and the solution from the medium-level adapted 
grid shows a peak energy flux 4% higher than its progenitor. Also, the gra- 
dients captured by the adapted grids are slightly steeper than their respec- 
tive original cases, and occur over a narrower band. The full-width/half- 
maximum value of <pue> for the coarse adapted grid is 86% of the original 
coarse-grid solution, and that for the medium adapted grid is 96% of the 
unadapted medium-resolution result. In both peak flux and flux width, the 
trend of data is toward the high-resolution result when going from original 
to adapted-grid solutions. 

Considerations for Future Work 

Although we have addressed the rules mentioned in the Discussion 
section, we have made no provision guaranteeing an adequate number of 
particles per cell for proper collision mechanics. Additionally, for optimum 
use of computational resources, we desire roughly equal numbers of parti- 
cles per cell. This rule conflicts with cell resolution guidelines, since the 
local cel) population is proportional to the product of local density and cell 
volume. The conflict can be resolved somewhat by using variable particle 
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scale factors, a concept already extended by Oly nick to the limit, where 
each cell is treated as a separate region. K It, may be desirable to combine 
current adaptive mesh procedures with such an approach. 

The last, example illustrates that in adapting to an initial solution with 
grossly inadequate resolution, the level of improvement, is constrained by 
fixing the number of cells and regions to those of the initial solution. It 
would be useful to automatically increase or decrease cell resolution by 
changing the number of cell rows or columns allotted to a given region 
based on the initial solution. It may even be advantageous to relax the 
constraints on the locations of certain outer boundaries (ahead of a bow 
shock into the freesteam, for example). 

It would also be advantageous to adapt rows and columns of cells 
simultaneously for enhanced ability to capture governing flowfield phenom- 
ena, within the constraints imposed by proper collision mechanics. Finally, 
restarting calculations from existing solutions on adapted grids, rather than 
beginning new computations from scratch, would greatly reduce the amount 
of time spent achieving new steady-state solutions. 

Conclusions 

We have begun developing a technique for adapting DSMC meshes to 
flowfield solutions. The functions used for adapting these grids are based 
on length scales that are physically meaningful in terms of satisfying certain 
empirical rules usually fulfilled manually. The formulation of these adap- 
tation functions is flexible enough that they may be used for unbounded 
flowfields as well as flows past bodies. Use of this technique can improve 
resolution without requiring greater resources in terms of numbers of cells 
and particles, or the extra CPU time and memory necessary to solve a 
larger simulation. 
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